% Function to run all model fit exercises. 
% Called from the Main script. 
% The six plots for Figure 1 are saved in the Output/ folder unless "outputfolder" is changed. 
% The function also outputs additional statistics that are referenced in the text, in matrix form (Modelfit)
function[Modelfit] = modelfit_func(PAR,Results,outputfolder) 
    
    % For reproducability when drawing random draws from vector (to use with randample(s,..) i.e. for Qhat)
    s = RandStream('mlfg6331_64');
    format short g; %to avoid trailing zeros

    % Also using fixed simulation draws in PAR.v0probs_fit (1000 draws, whereas PAR.v0probs
    % has fewer draws, used in estimation)
    % And using PAR.values_gridpoints, also to fix the randomness
    PAR.sharenobidders = Results.sharenobidders;
    
    % Load Estimation data
    if PAR.highvalue == 0
        data      = readtable('Estdata.csv','TreatAsEmpty','NA');
        data_nor  = data(data.hasreserve==0,:);
        data_posr = readtable('Estdata_posr.csv','TreatAsEmpty','NA');
    else
        data      = readtable('Estdata_high.csv','TreatAsEmpty','NA');
        data_nor  = data(data.hasreserve==0,:);
        data_posr = readtable('Estdata_posr_high.csv','TreatAsEmpty','NA');
    end
    
    % Estimated bidder and seller value distributions
    %Generalized gaussian
    thetab=[Results.mub,Results.kb,Results.shapeb];
    thetas=[Results.mus,Results.ks,Results.shapes];
    CDFb = @(x) ggd_cdf(x,thetab);
    PDFb = @(x) ggd_pdf(x,thetab);
    CDFs = @(x) ggd_cdf(x,thetas);
    PDFs = @(x) ggd_pdf(x,thetas);

    PAR.maxv = fminbnd(@(x) (CDFb(x)-0.999)^2,thetab(1)-3,thetab(1)+3);
    PAR.minv = fminbnd(@(x) (CDFb(x)-0.001)^2,thetab(1)-3,thetab(1)+3);
      
    location=thetab(1);
    scale=thetab(2);
    shape=thetab(3);
    if shape > 0
         maxv = location+(scale/shape)-0.001;
         PAR.maxv=min(PAR.maxv,maxv);
    else %shape<0
         minv = location+(scale/shape)+0.001; 
         PAR.minv=max(PAR.minv,minv);
    end

   % if thetas(3)>0
   %     lb = -3*thetas(1);
   %     ub = min(10,thetas(1)+thetas(2)/thetas(3));
   % else
   %     lb = max(-10,thetas(1)+thetas(2)/thetas(3));
   %     ub = 3*thetas(1);
   % end
   

    %%%%%%%% MODELFIT PLOTS 

    % --- Observed vs predicted conditional bidder values  --- %
    % Compute empirical CDF of idiosyncratic bidder value, NPFvSH, on grid, NPvgrid
    [NPFvSH, NPvgrid] = ecdf(data_nor.Uhat);
    ugrid             = linspace(min(NPvgrid),max(NPvgrid),PAR.values_gridpoints);
    
    % Calculate cdf sec-highest value conditional on different number of bidders
    NPFvn = ones(PAR.values_gridpoints,length(unique(data_nor.numberbidders))-2);
    for i = 1:length(unique(data_nor.numberbidders))-2
        j = 2+i-1;
        NPFvn(:,i) = my_ecdf(data_nor.Uhat(data_nor.numberbidders==j),ugrid);
    end
    
    % Translate into nonparametric cdf values for different number of bidders (inversion Athey & Haile 2002)
    if PAR.highvalue==0
        lown = 1;
        highn = 10;
    else
        lown = 2;
        highn = 9;
    end
    FvNP = zeros(PAR.values_gridpoints,highn-lown+1);
    for j = lown:highn
        n = j+1;
        for i = 1:PAR.values_gridpoints %loop over ugrid
            FvNP(i,j-lown+1) = fzero(@(Fv) NPFvn(i,j) - (n*Fv^(n-1) - (n-1)*Fv^n), [0,1]) ;
        end
    end
    
    % Plot for different number of bidders
    mar = ['+','o','*','x','s','^','p','h','<','>','d','v'];
    figure
    hold on
    for i=1:highn-lown+1
        plot(ugrid,FvNP(:,i),'LineWidth',1,'color','#767676','Marker',mar(i),'DisplayName',sprintf('n=%d',2+i-1))
    end
    plot(ugrid,CDFb(ugrid),'LineWidth',5,'color','black','LineStyle','--',...
        'DisplayName',sprintf('Predicted'))
    ax = gca;
    ax.FontSize = 12;
    xlabel('Log values (homogenized)')
    ylabel('CDF')
    if PAR.highvalue==1
            legend('show','Location','southeast')
    else
    legend('show','Location','northwest')
    end
    axis square
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    lgd.FontSize = 14;
    if PAR.highvalue==0
        print([outputfolder,'Figure1a'],'-dpng')
    end


    % --- Observed vs predicted conditional seller values --- %
    % For positive reserve price auctions
    % That are used in estimation:
    % - Not a negative mark-up
    % - Estimated quality trimmed with PAR.trimquality
    % Homogenized
    
    % Data used for fitting
    u0_idx     = data_posr.U0hat>=Results.screeningvalue&data_posr.U0hat<=Results.u0star;
    u0grid     = linspace(min(data_posr.U0hat(u0_idx)),max(data_posr.U0hat(u0_idx)),PAR.values_gridpoints);
    NPFv0      = my_ecdf(data_posr.U0hat(u0_idx),u0grid);
    Fv0_estim  = (CDFs(u0grid)-CDFs(Results.screeningvalue))./(CDFs(Results.u0star)-CDFs(Results.screeningvalue));
    
    % Plot
    figure
    hold on
    plot(u0grid,NPFv0,'LineWidth',1,'color','#767676','DisplayName','Observed')
    plot(u0grid,Fv0_estim,'LineWidth',5,'color','black','LineStyle','--',...
        'DisplayName',sprintf('Predicted'))
    ax = gca;
    ax.FontSize = 12;
    xlabel('Log values (homogenized)')
    ylabel('CDF')
    legend('show','Location','NorthWest')
    axis square
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    lgd.FontSize = 14;
    if PAR.highvalue==0
        print([outputfolder,'Figure1b'],'-dpng')
    end



    % --- Observed vs predicted reserve prices --- %
    % For positive reserve price auctions
    % That are used in estimation:
    % - Not a negative mark-up
    % - Estimated quality trimmed with PAR.trimquality, u0hat trimmed with PAR.trimu0hat
    % Include expectation over quality ; non-homogenized
    R      = length(PAR.vprobs_fit);  % nr sims v (separately for 15 n, so less then v0probs_fit)
    Rv0    = length(PAR.v0probs_fit); % nr sims v0
    
    % Simulate values from estimatedated seller value distribution and from estimated quality
    U0draws  = arrayfun(@(i) fminbnd(@(x) ((CDFs(x)-CDFs(Results.screeningvalue))./(CDFs(Results.u0star)-CDFs(Results.screeningvalue))-PAR.v0probs_fit(i))^2,Results.screeningvalue,Results.u0star),1:Rv0)'; % Large sample of seller conditional values
    Qdraws   = randsample(s,data_posr.Quality(u0_idx),Rv0,'true');  % Large sample of product quality
    
    % Compute predicted reserve prices from V0draws and Qdraws
    rfunc = @(r,u0,Q) ((r- exp(u0+Q)/(1-Results.cS) - (1-CDFb(log(r)-Q))/(PDFb(log(r)-Q)) ).^2);
    rpred = arrayfun(@(u0,Q) fminbnd(@(r) rfunc(r,u0,Q),0,9999),U0draws,Qdraws);
    
    % Observed reserve prices / Quality
    robs = data_posr.reserve_bot(u0_idx);
    %Qobs = data_posr.Quality(u0_idx);
    %nobs = data_posr.numberbidders(u0_idx);
    
    % Plot
    bwidth = (max(robs)-min(robs))/50;
    figure
    hold on
    histogram(robs,'BinWidth',bwidth,'Normalization','pdf','DisplayName','Observed','EdgeAlpha',0,'FaceColor',[0.5 0.5 0.5])                                                            % Includes trimming Quantile at the (1-trim_obs)th percentile and excluding negative markups, done in estimation_func()
    histogram(rpred,'BinWidth',bwidth,'Normalization','pdf','DisplayName','Predicted','FaceAlpha',0,'EdgeAlpha',1,'EdgeColor','black','LineWidth',2)               % Predicted reserve prices including sample selection above: trimming Quantile and negative markups
    ax = gca;
    ax.FontSize = 12;
    xlabel('Values (non-homogenized)')
    ylabel('PDF ')
    legend('show')
    axis square
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    lgd.FontSize = 14;
    if PAR.highvalue==1
        xlim([0 600])
    else
        xlim([0 100]);
        print([outputfolder,'Figure1d'],'-dpng')
    end
    


    % --- Actual number bidders in r>0 auctions --- %
    % Simulate reserve (rpred from above), given Quality from data, N from Poisson (r>0), Tastes from Fv|Z
    % Draw number of auctions from simulated reserve (R<Rv0)
    % Same for Qpred
    rpred  = randsample(s,rpred,R);
    Qdraws = randsample(s,Qdraws,R);
    
    % Simulate pot bidder values
    maxn   = PAR.maxn;
    Usim   = zeros(R,maxn);
    for n = 1:maxn
        quantilesim = PAR.vprobs_fit(:,n); % use pre-simulated quantiles for consistency
        Usim(:,n)  = arrayfun(@(i) fminbnd(@(x) (CDFb(x)-quantilesim(i))^2,-10,10),1:R)'; % Large sample of conditional bidder values
    end
    
    % Simulate number bidders in positive reserve auctions
    Nsim = random('Poisson',Results.lambdastarposr,R,1);
    
    % Add quality draw to U
    Vsim   = exp(Usim + Qdraws);
    
    % In each simulated auction i, for all bidders j > Nsim(i), set U value to minimum.
    barv = min(min(min(Vsim),min(rpred)))-1;
    for i = 1:R
        maxni = Nsim(i);
        Vsim(i,maxni+1:end) = barv;
    end
    
    % Get predicted sh values and hammer prices
    Hval    = max(Vsim,[],2);
    Vsim(Vsim==Hval)=min(min(Vsim));
    SHval   = max(Vsim,[],2);
    hampred = zeros(R,1);
    hampred(Nsim>1)  = max(SHval(Nsim>1),rpred(Nsim>1));
    hampred(Nsim==1) = min(rpred(Nsim==1),Hval(Nsim==1));
    hampred(Nsim==0) = rpred(Nsim==0);
    
    % Nr simulations predicted outside samples:
    % Actually, sample selected on the basis of hammer price not per bottle
    % But estimates based on these results
    % Why highestbid_bot is max 200? If total for multiple bottles is max
    % 200, then per one bottle also max 200
    if PAR.highvalue == 0
        length(hampred(hampred>200))/R
        hampred = hampred(hampred<=200);
        %Asim = Asim(hampred<=200);
    else
        %    length(hampred(hampred<200|hampred>1500))/R
        %    hampred = hampred(hampred>200&hampred<=1500);
        %    Asim = Asim(hampred>200&hampred<=1500);
        %    data_posr.highestbid_bot = data_posr.highestbid_bot(data_posr.highestbid_bot>200&data_posr.highestbid_bot<=1500);
    end
    
    % Plot
    figure
    hold on
    histogram(data_posr.highestbid_bot,'Normalization','pdf','EdgeAlpha',0,'FaceColor',[0.5 0.5 0.5])
    histogram(hampred,'Normalization','pdf','FaceAlpha',0,'EdgeAlpha',1,'EdgeColor','black','LineWidth',3)
    legend('show',{'Observed','Predicted'},'Location','northeast')
    ax = gca;
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    axis square
    xlabel('Values (non-homogenized)')
    ylabel('PMF');lgd.FontSize = 14;
    if PAR.highvalue==0
        print([outputfolder,'Figure1e'],'-dpng')
    end

    % --- Poisson fit r=0
    % Sample from estimated distributions
    bins        = 0:max(data_nor.numberbidders)-2;
    estimatedprobs = pdf('pois',bins,Results.lambdastarnor);
    bidderspred  = randsample(s,0:max(data_nor.numberbidders)-2,100000,'true',estimatedprobs);
    
    % Plot
    figure
    hold on
    histogram(data_nor.numberbidders(data_nor.numberbidders<=max(data_nor.numberbidders)-2),'Normalization','pdf','EdgeAlpha',0,'FaceColor',[0.5 0.5 0.5])
    histogram(bidderspred,'Normalization','pdf','FaceAlpha',0,'EdgeAlpha',1,'EdgeColor','black','LineWidth',3)
    legend('show',{'Observed','Predicted'},'Location','northeast')
    ax = gca;
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    axis square
    xlim([-1 max(bins)+1])
    xlabel('Number bidders per listing')
    ylabel('PMF');lgd.FontSize = 14;
    if PAR.highvalue==0
        print([outputfolder,'Figure1f'],'-dpng')
    end


    % --- observed versus predicted bids. Including estimated Quality
    % Estimation sample :
    % no reserve price auctions and more than 1 bidder
    
    % Simulate values
    maxn   = max(data_nor.numberbidders);
    R      = length(PAR.vprobs_fit); %nr sims
    Udraw  = zeros(R,maxn-1);
    for n = 2:maxn
        quantilesim   = PAR.vprobs_fit(:,n-1); % use pre-simulated quantiles for consistency
        Udraw(:,n-1)  = arrayfun(@(i) fminbnd(@(x) (CDFb(x)-quantilesim(i))^2,-10,10),1:R)'; % Large sample of conditional bidder values
    end
    
    % Simulate number bidders | exceeding 1 in no reserve auctions
    Ndraw = randsample(data_nor.numberbidders(data_nor.numberbidders>1),R,1);
    
    % In each simulated auction i, for all bidders j > Nsim(i), set value to minimum.
    for i = 1:R
        maxni = Ndraw(i);
        Udraw(i,maxni+1:end) = min(min(Udraw));
    end
    
    % Select the second-highest value in each auction, add draw of Quality, and put on normal scale
    Vdraw   = Udraw;
    Hval    = max(Vdraw,[],2);
    Vdraw(Vdraw==Hval)=min(min(Vdraw));
    SHval   = max(Vdraw,[],2);
    Qdraw   = randsample(data_nor.Quality(data_nor.numberbidders>1),R,1);
    bidpred = exp(SHval+Qdraw);
    
    % Observed bids
    bidobs  = data_nor.highestbid_bot(data_nor.numberbidders>1);
    nobs_nor= data_nor.numberbidders(data_nor.numberbidders>1);
        
    % Plot - homogenized bids
    bidobs_hom  = log(data_nor.highestbid_bot(data_nor.numberbidders>1))-data_nor.Quality(data_nor.numberbidders>1); %Ush
    bidpred_hom = SHval;
    bwidth = (max(bidobs_hom)-min(bidobs_hom))/50;
    figure
    hold on
    histogram(bidobs_hom,'BinWidth',bwidth,'Normalization','pdf','DisplayName','Observed','EdgeAlpha',0,'FaceColor',[0.5 0.5 0.5])                                                            % Includes trimming Quantile at the (1-trim_obs)th percentile and excluding negative markups, done in estimation_func()
    histogram(bidpred_hom,'BinWidth',bwidth,'Normalization','pdf','DisplayName','Predicted','FaceAlpha',0,'EdgeAlpha',1,'EdgeColor','black','LineWidth',2)               % Predicted reserve prices including sample selection above: trimming Quantile and negative markups
    ax = gca;
    ax.FontSize = 12;
    xlabel('Values (homogenized)')
    ylabel('PDF ')
    legend('show')
    axis square
    ax.FontSize=14; ax.TitleFontSizeMultiplier = 1.2;
    lgd.FontSize = 14;
    if PAR.highvalue == 1
        xlim([4.8 6.8])
    else
        xlim([1 6])
        print([outputfolder,'Figure1c'],'-dpng')
    end
    

   %%%%%%%% END MODELFIT PLOTS 




    % --- Table with stats referred to in results section --- %
    % Absolute mean squared deviations highest bids by N
    Modelfit    = zeros(18,1);
    
    % - Absolute deviations highest bids as share of average highest bid
    for n = 2:10
        Modelfit(n-1) = (abs(mean(bidpred(Ndraw==n))-mean(bidobs(nobs_nor==n))))./mean(bidobs(nobs_nor==n));
    end
    Modelfit(10)  = abs(mean(bidpred)-mean(bidobs))./mean(bidobs);  % Expectation over n : bidpred and bidobs are already N-weighted

    % Observed vs Predicted reserves 
    [h_res,p_res,ksstat_res] = kstest2(robs,rpred);
    Modelfit(11) = p_res;

    % --- Poisson fit r=0 --- %
    % For auctions without a reserve price (where number of bidders is observed)
    % create bins, observed counts per bin
    bins        = 1:max(data_nor.numberbidders)-3;%range = to make sure don't have bins with <5 expected counts, which will be grouped by chi2gof
    estimatedprobs = pdf('pois',bins,Results.lambdastarnor);
    obscounts   = arrayfun(@(x) length(data_nor.numberbidders(data_nor.numberbidders==x)), bins);
    expcounts   = sum(obscounts).* estimatedprobs;
    
    % Chi-squared test with as Null that observed counts generated from
    % Poisson, h = {0,1} = {cannot reject at 5% level, can reject at 5% level}
    [h,p_pois,st]    = chi2gof(bins,'Ctrs',bins,'Frequency',obscounts,'Expected',expcounts,'NParams',1);
   
    % KS test Poisson distribution fit with predicted (no reserve price)
    Modelfit(12) = p_pois;
   
    % To describe unscaled listing inspection costs
    %Estimating quality without N dummies
    % r>0
    Ndums_posr    = [data_posr.N3,data_posr.N4,data_posr.N5,data_posr.N6,data_posr.N7,data_posr.N8,data_posr.N9,data_posr.N10,data_posr.N11,data_posr.N12,data_posr.N13];
    Ncoefs        = [Results.LMcoef(end-10:end)];
    Nres_posr      = Ndums_posr*Ncoefs;
    data_posr.QnoN = data_posr.Quality-Nres_posr; 
  
    % Check unscaled listing inspection costs as share hammer price per bottle 
    unscaled_posr=(exp(data_posr.QnoN).*Results.eBo_posr)./data_posr.highestbid_bot;
    unscaled_posr=unscaled_posr(data_posr.numberbidders>1);
  
    % Record Median unscaled listing inspection cost r>0
    Modelfit(13) = quantile(unscaled_posr,0.5);

    % Same but r=0
    Ndums_nor     = [data_nor.N3,data_nor.N4,data_nor.N5,data_nor.N6,data_nor.N7,data_nor.N8,data_nor.N9,data_nor.N10,data_nor.N11,data_nor.N12,data_nor.N13];
    Ncoefs        = [Results.LMcoef(end-10:end)];
    Nres_nor      = Ndums_nor*Ncoefs;
    data_nor.QnoN = data_nor.Quality-Nres_nor; 
    unscaled_nor=(exp(data_nor.QnoN).*Results.eBo_nor)./data_nor.highestbid_bot;
    unscaled_nor=unscaled_nor(data_nor.numberbidders>1);
    
    % Record Median unscaled listing inspection cost r=0
    Modelfit(14) = quantile(unscaled_nor,0.5);

    % Mean bidder and seller values in population and on platform
    mean_V   = Results.mub - (Results.kb/Results.shapeb)*(exp((Results.shapeb.^2)/2)-1);
    mean_V0  = Results.mus - (Results.ks/Results.shapes)*(exp((Results.shapes.^2)/2)-1);
    Modelfit(15)=exp(mean_V);
    Modelfit(16)=exp(mean_V0);
    
    % Simulate values for sellers on the platform
    U0draws         = arrayfun(@(i) fminbnd(@(x) ((CDFs(x)-CDFs(Results.screeningvalue))./(CDFs(Results.u0star)-CDFs(Results.screeningvalue))-PAR.v0probs_fit(i))^2,Results.screeningvalue,Results.u0star),1:length(PAR.v0probs_fit))'; % Large sample of seller conditional values
    Modelfit(17)=exp(mean(U0draws));
   
    % Some other statistics in estimation sample for in the text;
    Modelfit(18) = mean(data_posr.highestbid_bot);

    % Save as a table
    Modelfit = table(Modelfit,'RowNames',...
        {'Winbid_n2','Winbid_n3','Winbid_n4','Winbid_n5','Winbid_n6','Winbid_n7','Winbid_n8','Winbid_n9','Winbid_n10',...
        'Winbid_expected','pvalue_KStest_reserve','pvalue_KStest_poisson',...
        'Unscaled_med_listing_rpos','Unscaled_med_listing_r0','MeanV','MeanV0','MeanV0platform','Avg_highestbid_bot_sample'});


end